function F = objfun105(x,feedback,settings)
% F = objfun105(x,feedback,settings)
% Cost function of MPC105 which controls mean(R^2) and var(R^2):
%    F =
%    sum_{i=1}^{i=P}(((R2_set-R2_i)/R2_set)^2)
%
% var(R^2) is from analytical solution
% mean(R^2) is calculate using analytical solution of alpha^2 and beta^2

% Identified parameters based on fitting of mean(R2) for 100s

model_R2_settings.W      = [0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50];
model_R2_settings.nu     = [4.788060e-3,3.107336e-3,1.243365e-3,3.408107e-4,2.668339e-4,1.687947e-4,1.322633e-4,7.475083e-5,6.448263e-5];
model_R2_settings.sigma2 = [1.708136e-2,1.688841e-2,1.277844e-2,8.336605e-3,8.988716e-3,8.874020e-3,9.026618e-3,8.661543e-3,9.375721e-3];     

if nargin == 2
    varR2_set   = 0;
    Fact_varR2  = 0.5;
    R2_set      = 50;
    Fact_r2     = 0.5;
    P           = 5;  % Prediction steps
    m           = 20; % Mode
    dt          = 10;
elseif nargin == 3
%     assert(isstruct(settings),'\nsettings should be a struct');
    varR2_set   = settings.varR2_set;
    Fact_varR2  = settings.Fact_varR2;
    R2_set      = settings.R2_set;
    Fact_r2     = settings.Fact_r2;
    P           = settings.P;
    m           = settings.m;
    dt          = settings.dt;
else
    error('\nThere must be two or three input');
end

% assert(isstruct(feedback),'\nfeedback should be a struct');
meanAlpha2 = (feedback.alpha).^2;
meanBeta2  = (feedback.beta).^2;
varAlpha2   = zeros(1,m);
varBeta2    = zeros(1,m);
meanR2      = zeros(P,1);   % Only predicted meanR2 is saved
varR2       = zeros(P,1);   % Only predicted varR2 is saved

for i = 1:P
%     if Tab_R_DepRate(1) <x(i) <= Tab_R_DepRate(end)
%         nu = interp1(Tab_R_DepRate,Tab_R_nu,x(i),'linear','extrap');
%         sigma2 = interp1(Tab_R_DepRate,Tab_R_sigma2,x(i),'linear','extrap');
%     elseif x(i) <= Tab_R_DepRate(1)
%         nu = Tab_R_nu(1);
%         sigma2 = Tab_R_sigma2(1);
%     else
%         nu = Tab_R_nu(end);
%         sigma2 = Tab_R_sigma2(end);
%     end
    
    % <R^2> model
%     for j = 1:m
%         temp = sigma2/(2*nu*j^2);
%         temp2 = exp(-2*nu*j^2*dt);
%         alpha2(i+1,j) = temp+(alpha2(i,j)-temp)*temp2;
%         beta2(i+1,j)  = temp+(beta2(i,j)-temp)*temp2;
%     end
%     meanR2(i) = sum(alpha2(i+1,:)+beta2(i+1,:))/(2*pi);

    % var(R^2) model
%     for j = 1:m
%         temp1 = exp(-2*nu*j^2*dt);
%         temp2 = exp(-4*nu*j^2*dt);
%         temp3 = sigma2/(nu*j^2);
%         temp4 = temp3*(2*temp2-2*temp1);
%         temp5 = 0.5*temp3^2*(temp2-2*temp1+1);
%         varAlpha2(j) = -alpha2(i,j)*temp4+temp5;  
%         varBeta2(j)  = -beta2(i,j)*temp4+temp5;
%     end
%     varR2(i) = sum(varAlpha2+varBeta2)/(4*pi^2);
    [meanAlpha2,meanBeta2,varAlpha2,varBeta2,meanR2(i),varR2(i)] = model_R2(x(i),meanAlpha2,meanBeta2,varAlpha2,varBeta2,dt,model_R2_settings);
end

F = 0;
for i = 1:P
    F = F+Fact_r2*((R2_set-meanR2(i))/R2_set)^2+Fact_varR2*(varR2(i))^2;
end